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Qh! We study the effects of including a running coupling constant in high-density QCD evo- 

' lution. For fixed coupling constant, QCD evolution preserves the initial dependence of the 

- ^ , saturation momentum Qg on the nuclear size A and results in an exponential dependence on 

^ ■ rapidity Y, Q'^{Y) = Q'^{Yq) exp [as d{Y — Yq)]. For the running coupling case, we re-derive 

analytical estimates for the A- and y-dependences of the saturation scale and test them nu- 
merically. The ^-dependence of Qs vanishes oc 1/VY for large A and Y. The ^-dependence 
is reduced to Q'^(Y) oc exp {A'\/Y + X) where we find numerically A' ~ 3.2. We study 
the behaviour of the gluon distribution at large transverse momentum, characterizing it by 
an anomalous dimension 1 — 7 which we define in a fixed region of small dipole sizes. In 
contrast to previous analytical work, we find a marked difference between the fixed coupling 
(7 ~ 0.65) and running coupling (7 ~ 0.85) results. Our numerical findings show that both 
a scaling function depending only on the variable r Qg and the perturbative double-leading- 
logarithmic expression, provide equally good descriptions of the numerical solutions for very 
small r-values below the so-called scaling window. 
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1 Introduction 



High-density QCD [1] - the regime of large gluon densities - provides an experimentally 
accessible testing ground for our understanding of QCD beyond standard perturbation 
theory. The Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation [2,3] is the perturbative 
framework in which the evolution of parton densities with decreasing Bjorken-x (in- 
creasing energy) is usually discussed. In the BFKL equation it is implicitly assumed 
that the system remains dilute throughout evolution and hence correlations between 
partons can be neglected. The fast growth of the gluon density predicted by the BFKL 
equation and experimentally observed at HERA, eventually leads to a situation in 
which individual partons necessarily overlap and, therefore, finite density effects need 
to be included in the evolution. These effects enter the evolution non-linearly, taming 
the growth of the gluon density. 

The need for and role played by saturation effects was first discussed in [4,5]. It 
was later argued [6-8] that in the high-density domain a hadronic object (hadron or 
nucleus) can be described in terms of an ensemble of classical gluon fields, and that 
the number of gluons with momenta smaller than the so-called saturation scale, is as 
high as it may be (i.e. saturated). The quantum evolution of the hadronic ensemble 
can be written in terms of a non-linear functional equation [9-15] where the density 
effects are treated non-perturbatively (see also [16,17]). 

An alternative approach, followed by Balitsky [18], relies on the operator product 
expansion for high-energy QCD to derive a hierarchy of coupled evolution equations 
(see [19] for a more compact derivation). In the limit of a large number of colours, the 
hierarchy reduces to one closed equation. This equation was derived independently by 
Kovchegov [20] in the dipole model of high-energy scattering [21-23]. 

The relation between these two approaches has been extensively discussed [13- 
15,24-27]. Apart from possible differences between the evolution equations in the 
kinematical region where the projectile becomes dense [24], the different approaches 
yield the same result, usually known as the Balitsky- Kovchegov (BK) equation. This 
equation has served as the starting point for a large number of analytical and numerical 
studies. It has also been derived in the ^'-matrix approach of [28] , and as the large- A^c 
limit of the sum of fan diagrams of BFKL ladders [29,30]. It corresponds, as BFKL, 
to a re-summation of the leading terms in a^ln (s/sq) (leading-log approximation). 
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Although the full analytical solution of the BK equation is not known, several of 
its general properties, such as the existence and form of limiting solutions, have been 
identified in both analytical [31-37] and numerical [29, 38-43] studies. Most of them 
refer to the fixed coupling case without impact parameter dependence, but analyses of 
the effect of a running coupling [42,44-48] and of the dependence on impact parameter 
[49-51] have also been carried out. Besides, there have been attempts to go beyond the 
large- Nc limit, either by analytical arguments [52-54] or by numerically solving the full 
hierarchy of evolution equations [47]. In this latter work, non- leading Nc corrections 
are found to give a contribution smaller than 10 -j- 15%, in qualitative agreement with 
what could be naively expected from a numerical correction of 0{1/N^). Prom a 
phenomenological point of view, studies of the BK equation are motivated by the 
geometrical scaling phenomenon observed in lepton-proton [55] and lepton-nucleus data 
[56, 57] which has been related to the scaling properties of the solution of the BK 
equation (see e.g. [58, 59] for recent numerical analyses of HERA data based on non- 
linear evolution). Further interest comes from the study of nuclear collisions [60], where 
saturation physics is argued [61] to underlie a large body of data including multiplicity 
distributions [57,62-66] and the rapidity dependence of the Cronin effect [43,67-70]. 

Next-to- leading- log contributions [71, 72] are known to have a strong impact on 
the BFKL equation [73-77]. Both the choice of scale in the coupling constant [78] 
and the implementation of kinematical cuts for gluon emission [42,79,80], together 
with physically motivated modifications of the kernel [81-83], have been proposed to 
mend some observed pathologies of next-to-leading-log BFKL. It is usually expected 
that the unitarity corrections included in the BK equation become of importance for 
parametrically smaller rapidities [74, 75] than those for which running coupling effects 
must be included [84]. This can only be definitively estabhshed once next-to-leading- 
log contributions are fully computed for BK (see [85] for a first step in this direction) . 
However, the inclusion of running coupling effects in BK may offer a hint of some of 
the effects induced at next-to-leading-log, as has been previously the case for BFKL. 
It may also help to reconcile the results of the equation with phenomenology [45,57]. 

In this paper we investigate numerically the influence of the running coupling on the 
solution of the BK equation without impact parameter dependence, leaving this last 
point for a future publication. We go beyond previous numerical studies [42,46,47] by 
making a detailed comparison between analytical estimates and our numerical solution 
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of the BK equation, and analyzing the Y- and A-dependence of the saturation scale. 
Our key results are the confirmation of the Y- and A-dependence of the saturation 
scale proposed analytically [32-34] , and the novel finding that the anomalous dimension 
(extracted for dipole sizes smaller than the inverse saturation scale) is different in the 
fixed and running coupling cases. To compare to analytical results which have been 
derived for asymptotically large energies, we shall evolve numerically to very large 
rapidities (up to F ~ 80), significantly beyond the experimentally accessible range. 

The plan of the paper is as follows. We first introduce the BK equation in Section 
121 and the different implementations of the running of the coupling constant in Section 
01 In Section 01 we explain the numerical method used to solve the BK equation. In 
Section El we present our numerical results, and we compare with previous numerical 
works and with analytical estimates. Finally, we summarize and discuss our main 
conclusions. 



2 The Balitsky-Kovchegov equation 

The BK equation gives the evolution with rapidity Y = ln(s/so) = ln(xo/x) of the 
scattering probability N{x,y,Y) of a qq dipole with a hadronic target, where x {y) is 
the position of the q (g) in transverse space with respect to the center of the target. 
We define 

- - - ? ^ + y ^ ^ ^ ^ ^ ^ 

r = X — y, = — - — , ri = X — z, r2 = y — z. [1) 

If one neglects the impact parameter dependence (which is justified for r <^ b, i.e. an 
homogeneous target with radius much larger than any dipole size to be considered), 
the BK equation reads (r = |r|) 

dN(r Y) r cPz 

^ ' - ' K{r, fi, fa) [iV(ri, Y) + Nir^, Y) - iV(r, Y) - N{r^, Y)N{r2, Y)\ , 



dY J 2n 

(2) 



where the BFKL kernel is 



K{f, fi, fa) = as , as = -^—^ ■ (3) 



1' 2 

The coupling constant is fixed and the kernel is conformally invariant. This implies 
that no impact parameter can be generated if not present in the initial condition. Also, 
there is no divergence for ri,r2 provided N{r,Y) oc for r with (3 > 0. 
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Figure 1: Diagrams for gluon emission in the evolution of a dipole and its Nc oo 
limit. 

This comes from the cancellation between real and virtual corrections inherited from 
the BFKL equation. The azimuthally symmetric form of the BFKL equation, which 
gives the dominant contribution at high energies, corresponds to Equation Q without 
the non-linear term. 

The BK equation has the following probabilistic interpretation [24] (see Figure^: 
when evolved in rapidity, the parent dipole with ends located at x and y emits a gluon, 
which corresponds in the large- A^^^ limit to two dipoles with ends {x, z) and (£, y) , 
respectively. The probability of such emission is given by the BFKL kernel (jS)), and 
weighted by the scattering probability of the new dipoles minus the scattering probabil- 
ity of the parent dipole (as the variation with rapidity of the latter is computed). The 
non-linear term is subtracted in order to avoid double counting. It is this non-linear 
term which prevents, in contrast to BFKL, the amplitude from growing boundlessly 
with rapidity. The BK equation ensures unitarity locally in transverse configuration 
space, |A^(r, y)| < 1. This is guaranteed since for A^(r, F ) = 1, the derivative with 
respect to F in (0) cannot be positive. 



3 Running coupling 

The BK equation ((21) was derived at leading order in asln(s/so) for a fixed coupling 
constant a^- An important part of the next-to- leading- log corrections is expected to 
come, as in BFKL, from the running of the coupling. The scale of the running coupling 
can only be determined when the next-to-leading-log calculation is available. In this 
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paper we introduce heuristically the running of the couphng, as done previously in 
BFKL (see e.g. [86,87]); we will use different prescriptions for the scales in order to check 
the sensitivity of the results. To motivate our choices, we recall the interpretation of 
the BFKL kernel as the Weizsacker- Williams probability for gluon emission written 
in a dipolar form. 



K{r,ri,r2) 



Ko{r,ri,r2) 



TT rrr. 



2^2 



1' 2 



4vr2 



gsn gsr2 



r 



(4) 



with Qs 

Three distance scales appear in this kernel: an 'external' one, the size of the parent 
dipole, r, and two 'internal' ones, the sizes of the two newly created dipoles, ri and r2- 
The latter depend on the transverse position of the emitted gluon z and on r through 
(P). We study three different prescriptions for implementing these scales in a running 
coupling constant in the BFKL kernel (jH): 



1. In the first modified kernel, Kl, the scale at which the running of the coupling 
is evaluated is taken to be that of the size of the parent dipole, r. This choice 
amounts to the substitution as in Equation (jH), 

as{r)Nc 



Kiir,fi,r2) 



TT 



(5) 



2. To implement the running of the coupling at the internal scale, we alternatively 
modify the emission amplitude in ^ before squaring it, 

2 



K2(f,ri,f2) 



47r2 



9s[ri)ri gs[r2)r2 



To 



(6) 



3. In order to check the sensitivity of the results to the Coulomb tails of the kernel, 
we further modify the kernel K2 by imposing short range interactions, so that 
the emission of large size dipoles is suppressed. To do this, we weight the gluon 
emission vertex by exponential (Yukawa-like) terms. 



K3{f,ri,r2) 
with /i = Aqcd- 



47r2 



-/ir2/2 



9s{r2)r2 



(7) 



Let us anticipate that the different prescriptions Kl, K2 and lead to very 
similar results for the evolution. This can be traced back to the fact that all the 
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geometrical dependence on z is integrated out so that only the r dependence in the 
running of the coupling survives. Even the introduction of the exponential damping 
has little effect, unless the range of the interaction is chosen unphysically small (i.e. 
\i ^ Kqcd)- However, the inclusion of a short range damping effect is known [49,50] to 
alter significantly the solution of the BK equation with impact parameter dependence, 
which we do not consider in the present work. 

For the qualitative properties of BK evolution studied in this paper, the precise 
value and running of the coupling constant is unimportant. To be specific, we use the 
standard one loop expression 



where A is an infrared regulator and /?o = ~ 2A^j with N j = 3. Both A and 

^QCD are determined from the conditions as{r = oo) = ao, as{r = 2/M^o) = 0.118, 
where M^o is the mass of the boson. In our work, this choice is not motivated 
by phenomenology, but by its use in related works e.g. [32, 45] to which we want to 
compare. From now on, when comparing fixed and running coupling results, it will be 
understood that the value for the fixed coupling is the same as the one at which the 
running coupling is frozen, Oq- 

4 Numerical method and initial conditions 

To solve the integro-differential equation we employ a second-order Runge-Kutta 
method with a step size AY = 0.1. We discretize the variable |r| into 1200 points 
equally separated in logarithmic space between r^m = 10^^^ and Vmax = 10^. The nu- 
merical values of these limits are dictated by the initial conditions and Aqcd- Through- 
out this paper, the units of r will be GeV^^, and those of Qg GeV. The integrals in 
Q are performed with the Simpson method. Inside the grid a linear interpolation is 
used. For points lying outside the grid with r < r^m a power-law extrapolation is used, 
while for points with r > Vmax the saturated value of the scattering probability is held 
constant, N{r) = N{rmax) = 1- While the initial conditions of N{r) give negligible 
values for r small but much larger than r^j^, the evolution leads to a gradual filling of 
values close to r^m with increasing rapidity, which would result eventually in numerical 
inaccuracies. To solve this problem and push the evolution to very large rapidity, we 




127r 
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rescale, in the fixed coupling case, the variable r in the solutions at intermediate values 
of Y and use them as initial condition (a power-law extrapolation is used for small 
values of r in order to cover the r-range lost in the rescaling procedure). In this way, 
we are able to evolve initial conditions with ~ 1 GeV up to F ~ 36 for = 0.4 
and up to y ~ 72 for = 0.2. In the running coupling case, the evolution is much 
slower and this rescaling is not needed to get to large rapidities. The accuracy of our 
numerical solution for all r-values inside the grid is better than 4% up to the largest 
rapidities. It is much better than 4% in most of the r-region studied. We have checked 
this numerical accuracy by varying the step size in Y, by comparing our results to those 
of a fourth-order Runge-Kutta method, by varying the limits of the grid, by doubling 
the number of points used to discretize the function in the grid, and by using different 
integration, extrapolation and interpolation methods. 

We evolve three different initial conditions starting from some fixed value of xq 
(in practice one usually takes Xq ~ 0.01). The first initial condition we refer to as 
GBW since it shows at fixed Xq the same r-dependence as the Golec-Biernat-Wiisthoff 
model [88]: 



N^^^ir) = 1-exp 



12 



(9) 



4 

However, in contrast to the GBW model [88], our x-dependence comes from BK evo- 
lution and we do not impose a power-law parameterization of the x-dependence of Q'^. 
Here and in the other initial conditions (fTUj) . (fTT|) below, we denote as Q'^ what is usu- 
ally called the saturation scale. Our definition of the saturation scale Qs is somewhat 
different, see Equation (fT^ below, but the relation between both scales is straightfor- 
ward e.g. in GBW, = —4 In (1 — k) Q^. The second initial condition takes the form 
given by the McLerran-Venugopalan model [6,7] (MV): 

r^Qf, / 1 



Ar^^(r) = 1 - exp 



4 '-'te^^ 



(10) 



These initial conditions have been used in previous works e.g. [39,43]. For transverse 
momenta ~ 1/r > 0{\ GeV), the sensitivity to the infrared cut off e is negligible. The 
amplitudes N'^^'^ and A^*^^ are similar for momenta of order Q'^ but differ strongly 
in their high-/c behaviour. The corresponding unintegrated gluon distribution (l){k) = 
J ^^e***'^ A^(r) decays exponentially for N^^^ but has a power-law tail ~ for 
]\jMV ^ As a third initial condition, we consider 

Ar^^(r) = l-exp[-(rQ'J'=]. (11) 



The interest in this ansatz is that the small-r behaviour A^^"^ oc corresponds to an 
anomalous dimension 1 — 7 = 1 — c/2of the unintegrated gluon distribution at large 
transverse momentum. This anomalous dimension can be chosen to differ significantly 
from that of the initial conditions N'^^'^ and N^^ . Our choices c = 1.17 and c = 0.84 
are somewhat arbitrary. They can be motivated a posteriori by the observation that the 
anomalous dimension of the evolved BK solution for both fixed and running coupling 
lies between the anomalous dimension of the initial conditions A^"^"^ and A^<^^^ (or 
N^^^). Thus, the choice of A^^"^ is very convenient to establish generic properties of 
the solution of the BK equation. The values of Q'^ in Equations 0, (fTUI) and pi|) are 
1.4 GeV for GBW, 4.6 GeV for MV, 0.7 GeV for AS with c = 1.17 and 0.6 GeV for 
AS with c = 0.84. These values have been used in all our studies except in those on 
the A-dependence in Section 15. 4t where Q'^ has been rescaled with the nuclear size as 
discussed in that Section. 

5 Results 

In this Section, we discuss our numerical results and how they compare to previous 
numerical work and analytical estimates. 

5.1 Evolution: Insensitivity to details of running coupling pre- 
scription 

Figure El shows the evolution of the dipole scattering probability for GBW initial con- 
dition with fixed and running coupling. The evolution is much faster for fixed coupling 
than for running coupling, as already known from previous numerical studies [42,46,47]. 
Remarkably, the solution is rather insensitive to the precise prescription with which 
running coupling effects are implemented in the modified BFKL kernels Kl, K2 and 
K?). These differences are very small compared to those between fixed and running 
coupling. 

The small differences arising from the use of different kernels can be understood 
qualitatively. For example, compared to Kl, the results obtained for K2 are enhanced 
at small values of r and suppressed at large values of r. This is due to the fact that 
e.g. for a typical size ~ l/Qs of the emitted dipoles ri, r2, a larger size r > l/Qs of 
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Figure 2: Solutions of the BK equation for GBW initial condition (dotted line) for 
rapidities Y = 6, 12 and 18 with ao = 0.4. Left plot: Evolution with fixed {KO, solid 
lines) and running coupling {Kl, dashed lines). Right plot: evolution with running 
coupling for kernel modifications Kl (solid lines), K2 (dashed lines) and K3 (dashed- 
dotted lines). 



the parent dipole amounts to a larger coupling gs{r) entering the kernel Kl than the 
couplings Qsi'^^i), gs{f2) entering K2. Thus, at large r the evolution is slower for K2, 
which results in the observed relative suppression. The analogous argument implies a 
relative enhancement obtained from the kernel K2 for small r < 1/Qs- 

Figure El also shows that the effects of imposing short range interactions, i^S, are 
very small (unless the range of the interaction is unphysically small). As expected, 
effects from short range interactions included in K?) are larger for larger values of r. 
It is conceivable that the main next-to-leading-log effects on the original BK kernel 
are those of the running of the coupling constant included here, and that further 
modifications, such us kinematical constraints [42,79,80], are comparatively small [89]. 
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5.2 Scaling 



In the limit Y ^ oo, the solutions of the BK evolution are no longer functions of the 
variables r and Y separately, but instead they depend on a single scaling variable 



Here, the saturation momentum Qs{Y) determines the transverse momentum below 
which the unintegrated gluon distribution is saturated. It can be characterized by the 
position of the falloff in N{r), e.g. via the definition 



where k is a constant which is smaller than, but of order, one. We have checked that 
different choices such as k = 1/2 and k = 1/e lead to negligible differences in the 
determination of Qs{Y). The results given below have been obtained for k = 1/2. 

In the fixed coupling case, the scaling property N{r, Y) — > N{t) has been quantified 
in previous numerical works [39,41,43] and confirmed by analytical calculations [35-37]. 
In the running coupling case, the scale invariance of the BFKL kernel is broken by the 
scale Aqcd and it is a priori unclear whether scaling persists. However, when the two 
scales in the problem are separated widely due to evolution to large rapidity, QsiX) ^ 
Aqcd, one may expect that the scaling property of the BK solution is restored. In 
agreement with previous numerical works [46,47], we confirm this expectation: for all 
modifications Kl, K2 and K3 of the BFKL kernel, the solutions tend to universal 
scaling forms as rapidity increases. Moreover, with increasing rapidity the sensitivity 
to the choice of scales in the kernel and its short range modification, as well as to the 
initial condition and to the value of the coupling constant in the infrared, becomes 
eventually negligible (see Figure E}. 

As seen in Figure El the shape of the scaling solution differs significantly for fixed 
and running coupling as observed already in [46]. The running of the coupling sup- 
presses the emission of dipoles of small transverse size (i.e. small r and large transverse 
momenta). This leads to an enhancement in the large r region of N{t) which is seen 
for the running coupling case in Figure 01 

The accuracy of scaling at small r has been studied in a previous work [43] for 
the fixed coupling case. Here we check scaling for both fixed and running coupling 



r = rQ,{Y). 



(12) 



iV(r = l/Q,(F),F) = «:, 



(13) 
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Figure 3: Scaling solutions of BK for y = 0, 20, 30 and 40 (plots on the left) and F = 0, 
40, 60 and 80 (plots on the right). Upper-left: evolution for fixed (solid) and running 
coupling (-ft'l, dashed lines) for GBW initial conditions. Upper-right: solutions for 
the kernels K\ (solid), KI (dashed) and K?) (dashed-dotted lines). Lower-left: scaling 
function for K\ with two different values of frozen coupling, ao = 0.4 (solid) and 
= 0.2 (dashed lines). Lower-right: scaling solutions with running coupling (-ft'l) for 
two different initial conditions, GBW (solid) and MV (dashed lines). In all plots the 
initial conditions correspond to the dotted hues and ao = 0.4 unless otherwise stated. 

by comparing our numerical results to the scaling forms proposed in Reference [33]. 
There, it was argued that in the so-called scaling window r^^ < r < 1, the asymptotic 
solution of A^(r, Y) takes the following scaling forms for fixed and running coupling, 
respectively [33]: 

/i)(r) = ar^^ (in r2 + 5) , (14) 

/2)(^)=ar2^|^lnr2 + i^ . (15) 

Here, 1 — 7 is usually called the anomalous dimension which governs the leading large- 
h behaviour of the unintegrated gluon distribution. We define 7 from a fit of our 
numerical results to the functions ()14j) and ()15|1 in the F-independent region 10^^ < 
r < 10~^, i.e. for 10^ > 1/t > 10 Qs, with a, 7 and 5 as free parameters. The 
results given below were found to be insensitive to a variation of the lower limit of this 
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Figure 4: The rapidity dependence of the parameter 7, characterizing the anomalous 
dimension 1 — 7, as determined by a fit of (|T4|) to the BK solutions for different initial 
conditions: GBW (squares), MV (circles), and AS with c = 1.17 (stars) and c = 0.84 
(triangles). Left plot: results for fixed coupling with Oq = 0.2. Right plot: results for 
running coupling with = 0.4 and two versions of the kernel Kl (empty symbols) 
and K2 (filled symbols). 

fitting range. 

For the case of fixed coupling constant, we find that the function /^^ provides a 
very good fit to the evolved solutions. In Figure |31 we show the fit values of the 
parameter 7, obtained for fixed coupling constant from the evolution of different initial 
conditions N'^^^ , A^*^^, and A^'^'^ for different values of c. At initial rapidity, these 
distributions have widely different anomalous dimensions but evolution drives them to 
a common value, 7 ~ 0.65, which lies close to the theoretically conjectured one [32,33] 
of 0.628. For a small fixed coupling constant = 0.2, this asymptotic behaviour is 
reached at F ~ 70, while for a larger coupling constant oq = 0.4 the approach to 
this asymptotic value takes half the length of evolution (results not shown). For fixed 
coupling solutions, /^^ does not provide a good fit to our numerical results. 

We have repeated this comparison for all running coupling solutions. We found that 
both and /^^ provide good fits and yield very similar values of 7. The results for 
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K3 are numerically indistinguishable from those for K2 and will not be shown in what 
follows. Also, the value of 7 was found to be independent of the coupling constant ao 
at r ^ 00. In Figure|31 we show the values of 7 extracted from a fit to Irrespective 
of the initial condition, they approach a common asymptotic value 7 ~ 0.85. While our 
numerical findings for A^^"^ with c = 0.84 are not inconsistent with the approach to this 
asymptotic value, no firm conclusions can be drawn. This initial condition just starts 
too far away from the asymptotic scaling solution to reach it within the numerically 
accessible rapidity range. In this case, the monotonic increase of 7 with rapidity at 
large Y is smaller than the increase for A^'^'^ with c = 1.17 at comparable values of 7, 
indicating that the rapidity evolution of the anomalous dimension depends in general 
not only on the small-r behaviour, but on the full shape of the scattering probability. 

The value 7 ~ 0.85 is considerably larger than the one found in fixed coupling 
evolution. This is in agreement with previous numerical results [46] but in contrast to 
theoretical expectations [32,33,45] which predict the same value of 7 for the fixed and 
running coupling cases. As an additional check, we have performed running coupling 
evolution from an initial condition given by the solution at large rapidity of fixed 
coupling evolution (for which 7 ~ 0.65). We find that even with this initial condition, 
running coupling evolution leads to a value of 7 ~ 0.85. 

It has been argued [32, 33] that expressions (fT^ and (fT3j) are only valid for values 
of T inside the scaling window, ~ Aqcd/Qs(X) < t < 1 with Yq the initial rapidity, 
and that the dipole scattering probability returns to the double-leading-log (DLL) 
expression 



Ar^^^(r) = a{Y) [- In (r^A^)]-^/^ exp \b{Y)^~\n {r^A^ 



(16) 



with a(Y) oc y^/^ and b(Y) oc VY, for values r < t^^. We have checked that this 
form provides a good fit (fit and numerical solution differ by less than ±10%) to the 
fixed coupling solution of BK for r < Ts„ = A/Qs(Y), A ~ 0.2 GeV, see Figure El 
Our comparison is limited to rapidities Y < 20, since the scaling window starts to 
extend over the entire numerically accessible r-space for Y > 20. Up to y = 20, 
the coefficients a{Y) and b{Y) follow the expected DLL F-behaviour, see Figure 
However, the scaling ansatz f^^ provides an equally good fit to the BK solutions for 
r < Tsw. This is the reason why in previous numerical studies [43] no upper bound for 
a scaling window was found. When the solutions of BK are fitted to f^^ within the 
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Figure 5: Plot on the left: solutions of the BK equation (solid lines) with GBW 
initial condition and fixed coupling = 0.2 compared to fits (dashed lines) to the 
DLL expression (|T6|l . for rapidities y = 0, 4, 8, 12, 16 and 20 (curves from right to 
left). Plots on the right: values of the coefficients a{Y) and b(Y) (circles) in the DLL 
expression versus Y, compared to fits (curves) to the functional form suggested by 



DLL. 

scaling window, the values of 7 at y = for both initial conditions are < 20% smaller 
than those found when the fit is done within a fixed r-window. But for larger Y the 
values of 7 extracted from fits within either the scaling window or some fixed r-window 
approach each other and quickly coincide. 



5.3 Rapidity dependence of the saturation scale 

In the scaling region, for large Y where QsiY) ^ ^qcd, the BK Equation ^ for 
fixed coupling constant can be written in terms of the rescaled variables r = Qs{Y)f, 
Ti = Qs(Y)fi and T2 = Qs(Y)r2- The F-dependence of N{r,Y) = N{t) is then 
contained in Qs{Y). Rewriting the derivative on the left-hand-side of (j21), 

dNjr) dQ,{Y) ON d\nmY)/A'] , ON 
dY dY ^ dr dY ^ dr^ ' ^ ^ 
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one finds [32] 

/ j^g^ ^ . i^(^) - ^(0)1 ^ . ^'^ . (18) 

Performing tlie same integration over cPr/r'^ = cPt/t'^ on tlie riglit-liand-side of 
one finds a number 

^ f^^^^^ + ^^^^^ " ^^""^ ~ Ar(r0iV(r2)], (19) 

wliicli is independent of Y . Tfie numerical value of d cannot be obtained witliout tlie 
knowledge of the scaling solution N{t), and several approximations have been proposed 
[32,33] which we will compare with our numerical results. Combining Equations 
p8|) and p9|) . the F-dependence of the saturation scale is determined [32] by 

dlnmY)/A'] 

gy = das. (20) 

Thus, for the case of a fixed coupling constant, the saturation scale grows exponentially 
with rapidity, 

Ql{Y) = Qlexp[AY], (21) 

where = = constant, A = da^ and Qq = Q^{Y = 0) (i.e. the evolution starts at 
Y = 0). 

For running coupling, the momentum scale is expected to be ~ Qs(Y). This sug- 
gests the substitution dg — > ds{Qs{y)) in Equation (^0]) . To see this explicitly, let us 
include, as in Kl, the coupling constant ds{r) in the integrand of (fT^ . which leads to 

12Nc fd^Td'^n 1 1 



Po J 2vr2 rf(r-fi)' In [Q^(>^)/A|cd]- In (rV4) 
X [N{n) + N{\r-n\) - N{r) - N{n)N{\r - (22) 

For r ^ 1, the integrand vanishes. For r ^ 1, the integral in d'^Ti is finite and 
the remaining d^r suppresses the contribution of small r. So we conclude that the 
dominant region is that of r ~ 1 and thus it is legitimate to approximate 

12Nc f d^rd^Ti 1 1 



/3o J 27r2 rl{f-f,y\nmY)/Klc^]-\r.{r^/A) 

X [N{T,)+N{\T-f,\)-N{T)-N{r,)N{\T-f,\)] ^ dds{Qs{Y)). (23) 

This approximation is also supported by numerical results [29,39,42] which show that 
in momentum space the typical transverse momentum of the gluons is ~ Qs- Due to 
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Figure 6: The rapidity dependence of the saturation momentum Ql for fixed as = 0.4 
(thick sohd), fixed = 0.2 (thin sohd), and running couphng with = 0.4 for kernels 
Kl (dashed), K2 (dashed-dotted) and K3 (dotted hues). For each group, hues from 
top to bottom in the rightmost side correspond to initial conditions AS with c = 1.17, 
MV and GBW. 

the similarities in the evolution shown previously, this should also hold for other imple- 
mentation of the scale of the coupling constant such as K2 and K3. The logarithmic 
dependence of the coupling constant on QsiY) in (f^ . combined with Equation (pUj) . 
leads to [32] 

^ (24) 



Ql{Y) = exp [A'VyTx\ , 

where (A')^ = 24A'crf//5o and X = (A'y^ln iQl/A"^). This estimate indicates that the 
rapidity dependence of the saturation scale is much weaker for running than for fixed 
coupling constant. 

Figure El shows the F-dependence of for several initial conditions and different 
choices of ao, calculated for all the kernels considered in this work. The rise of Qs is 
much faster for fixed than for running coupling, as already observed in [42,44-48,83,89]. 

For fixed coupling constant, exhibits with good accuracy an exponential be- 
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haviour for high-enough values of Y. The value of the slope extracted from a fit to the 
function (PTj) is A ~ 1.83 for do = 0.4. As expected, for do = 0.2 this value is reduced 
by a factor two, A ^ 0.91. For the constant (jl9j) . we find d ~ 4.57, in agreement with 
previous numerical studies at very high rapidities [43] but slightly smaller than the 
theoretical expectation d = 4.88 [32,33]. In previous numerical studies [39,40,42], an 
even smaller value of ~ 4.1 was obtained. We have checked that this is due to the 
fact that the rapidity region for the fit in our case corresponds to much larger Y. 

For the case of a running coupling constant, an exponential fit can be done only for 
a very limited F-region. For example, for F ~ 10 we find a logarithmic slope ~ 0.28 for 
GBW or MV initial conditions with Qo ~ 1 GeV, in agreement with the results of [45] 
but smaller than the values found in [83] (see also [48,89]). The exponential function 
12111 is unable to fit the full Y-range. In contrast, the weaker rapidity dependence of 
fl24|l does provide a good fit in the full Y-iange. The fit to ()24j) yields A' ~ 3.2, while 
the theoretical expectation [32, 33] is slightly larger. A' = 3.6. We finally note that 
in [47] the F-derivative of In Q^Y) has been found numerically to be proportional to 



y <^s{Qs{y)) in a much more restricted range of Y. We have been unable to fit our 
results over the full y-range to the corresponding Y-dependence, Ql(Y) oc expF^/^. 

We have found very little sensitivity of the values of A and A' to the fitting region, 
provided Y was chosen large enough. Our fits typically started at F ~ 15 where the 
asymptotic behaviour is approached, and explored the highest rapidities numerically 
accessible. Also, our results for running coupling do not depend on the choice of the 
kernel Kl, K2 or K3, on the initial condition or on the value of do. However, the AS 
initial condition with c = 0.84 is not included in our study since it does not approach 
the asymptotics within the numerically accessible rapidity range. 

In References [33,35-37,45] sub-leading terms in the F-behaviour of Qs have been 
presented. A form of the type dliaQl{Y)/dY = dga - hY~^ + cY^'^/"^ / {2^/d~s) has 
been proposed in the fixed coupling case, with a = 4.88, b = 2.39 and c = 2.74. This 
function contains all terms for the F-evolution of the saturation scale that are universal 
i.e. independent of the initial condition (see also [90] for a comparison of solutions of 
BK to this functional form). The constant term corresponds to Equation (j2ip . We have 
used this functional form to fit the results of fixed coupling evolution on d In Ql{Y) / dY 
for different rapidity regions within Y = 5 -r AO (72) for dg = 0.4 (0.2), for the GBW 
and MV initial conditions. First, we have used our definition of the saturation scale 
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fjl3|) with K = 1/2. From a simple comparison to the proposed expression (using the 
theoretical coefficients provided in [37]), we are able to clearly identify in our numerical 
results the presence of the first two terms. On the contrary, the presence of the third 
term is disfavoured. Fitting our numerical results to the first plus second terms, the 
value of a we find, a ~ 4.9, is quite stable with respect to variations of the fitting region. 
It is higher than the value of d we extract with only the linear term ()21|) . d ~ 4.57, and 
closer to the theoretical expectation d = 4.88 [32,33]. In this two-parameter fit we get 
a value of 6 ^ 2.4-^2.5, varying slightly with the y-region of the fit. This value is quite 
close to the theoretical expectation 2.39. On the other hand, in a three-parameter fit 
the values of b and c we extract are very unstable (even changing signs) with respect to 
variations of the lower limit of the fitting region between Y = 5 and 20. We have also 
tried to get the value of c from a fit to d/dY[Y d\nQl{Y)/dY] = aga — cY^^^'^/ {4:^/a^). 
While we find again a value of a ~ 4.9, the value of c turns out to depend, as in the 
previous analysis, considerably on the fitted y-region. Secondly, we have used the 
definition of the saturation scale (fT!?|l but now with k = 0.01 (i.e. we define Qs in a 
point in which the dipole scattering probability is far from its unitarity limit). In this 
case, a simple comparison to the proposed expression using the theoretical coefficients 
provided in [37] allows to clearly identify in our numerical results the presence of the 
three terms. Still, a three-parameter ffi to our numerical results does not provide values 
of b and c stable with respect to changes in the ffiting region. This influence of the 
definition of the saturation scale on the determination of the sub-leading corrections 
to its y-behaviour is consistent with the finding in [90]. 



5.4 Nuclear size dependence of the saturation scale 

The nuclear size enters the initial condition. The question is whether the BK evolu- 
tion modifies or preserves this initial A-dependence. For realistic nuclei, the impact 
parameter is likely to have an important effect on this A- dependence. This has been 
examined partially in [50,51]. However, the question is already of interest for the case 
without impact parameter dependence [34,39,40], which we study here. 

Let us ffist assume some arbitrary A-dependence which we include in the initial 
condition by the rescaling factor 

— > hr'^ (25) 
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(this is true for GBW and AS initial conditions but not for MV due to the presence 
of the logarithm; however, the numerical results for the A-dependence obtained with 
MV initial conditions are, for all purposes, equivalent to those with GBW). Here, h 
contains the information about the nuclear size, and Equation (pUjl reads 

/ A 21 



(26) 



In the case of a fixed coupling constant, the dilatation invariance of the BK equation 
(PI) allows to scale out any nuclear dependence included in the initial condition. Thus, 
the A-dependence of the saturation scale is unaffected by evolution. To explore the 
case of a running coupling constant, we use the one-loop expression for and write 



Ql{y) 



A2 



exp 



(A')2y + In^ 



A2 



(27) 



Multiplying by h for the nucleus to undo the rescaling, and setting h = 1 for the proton, 
we get 



exp 



(A')2y + In^ 



If we assume the hierarchy 

(A')'F > In' 

so y4 ^ 1, we find 

In 



lhQliY = 0)] 




A2 





(A')2r + w 



Q'sjY = 0) 
A2 



hQl{Y = Oj 
A2 



> In' 



Ql{Y = 0) 
A2 



A2 



:28i 



(29) 



(30) 



Here, hQliY = 0) is the initial saturation momentum for the nucleus, and Equation 
dnOI) coincides with Equation (44) of [34] with (A')' as defined below Equation (see 
also [47,64,91] for related discussions). This result suggests that any information about 
the initial A-dependence of the saturation scale is gradually lost during evolution: albeit 
at extremely large rapidities, all hadronic targets look the same. Usually, one assumes 
an A^/^-dependence of the saturation scale for the initial condition [6,7], h oc A^^^. 
However, other ^-dependencies have been proposed e.g. [92]. 

Figure [7| shows that fixed coupling evolution preserves the A-dependence of the 
saturation scale irrespective of whether this dependence is oc A^^^ as for the GBW or 
MV initial conditions (which produces numerical results for the A-dependence which 
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Figure 7: Upper plot: QIa/Q% versus A^^'^ for initial conditions GBW {QIaO^ = 0) oc 
solid) and AS with c = 1.17 {QIa^X = 0) oc A^/^^, dashed lines); thick lines 
are the results for F = in the running coupling case and for all rapidities in fixed 
coupling; for running coupling, different rapidities F = 10, 40 and 72 (thin lines) are 
shown from top to bottom for each initial condition. Lower plot: QIaI {^^^^Qlp) versus 
Y for GBW with A^/^ = 3, 6, 10, 20 and 50 with the same line convention as the upper 
plot (the results for fixed coupling have been obtained for F < 36 and extrapolated 
as a constant equal to 1). In all plots uq = 0.4 and in the running coupling case the 
kernel Kl has been used. 



are very close to those obtained for GBW), or whether it differs from oc A^^^ due to 
an anomalous dimension included e.g. in the AS initial condition. On the other hand, 
running coupling evolution is seen to reduce the A-dependence with increasing rapidity. 
We find that if fitted in a wide rapidity range, the dependence of In [Q'^a(X) / Q%0^)] 
on y is ~ F""'^. However, for large values of A and Y, the decrease with increasing Y 
is oc 1 / and thus well described by (jHUj) [34] . 

Combining the rescaling argument based in with the observation that the DLL 
solution is approached for small r or large transverse momentum k, one is led to an 
interesting implication for the large-fc behaviour of the ratios of gluon densities in 
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nuclei over nucleoli (or central over peripheral nucleus) [43,68-70]. In fixed coupling 
evolution the rescaling of the initial condition trivially implies the same rescaling 
in the evolved solution, which we will consider to be DLL for sufficiently large k. Thus 
one gets for the ratio R of the gluon densities in transverse momentum space for nuclei 
over nucleon 

-3/4 

exp I b{Y) 



R 



In h 



In /A2) - In /i - J\n (A;^ /A2) 



(31) 



ln(A;VA2) 

This ratio tends very slowly to 1 for — > oo. We have checked that the results of this 
formula agree with the numerical computations in [43] and thus it provides justification 
to the apparent absence of a return to the coUinear limit, R = 1 aX k ^ oo, found in 
this reference for the largest studied /c-values. 



6 Conclusions 

The inclusion of a running coupling constant may be expected to account for impor- 
tant next-to-leading-log effects in the BK equation, as has been previously the case for 
BFKL. This motivates the present numerical study of the BK equation without impact 
parameter dependence. Our main results are insensitive to details of the implementa- 
tion of running coupling effects, the infrared regulation of the coupling constant and 
the choice of initial conditions which are evolved. They can be summarized as follows: 

1. The rapidity dependence of the saturation momentum is much faster for fixed 
coupling constant than for the running one, as observed previously [42,44,46,47]. 
It is well described by Q^iY) oc exp {agdY) for fixed coupling and by Q^iY) oc 
exp {A'^/Y + X) for running coupling. For large rapidities, we find d ~ 4.57 
which is slightly smaller than the theoretical expectation d = 4.88 [32,33]. For 
running coupling, we find A' ~ 3.2, slightly smaller than the expected value 
A' ^ 3.6 [32,33,45]. For a very limited region of Y, a fit to the exponential form 
Ql{Y) oc exp (DY) works even for running coupling, but it cannot account for 
the entire F-range. For the fixed coupling case, we have checked the existence 
of the sub-leading terms in the ^-dependence of the saturation scale proposed 
in [33,35-37]. As found in [90], their precise determination depends on the 
definition of the saturation scale. 
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2. For sufficiently large rapidity, the solution of the BK equation with fixed coupling 
is known to show scaling [39,41,43]. We confirm scaling for the running coupling 
case in agreement with [46,47]. The approach to the scaling solution is faster 
with fixed than with running coupling. 

3. As observed previously [46] and at variance with analytical estimates [32,33,45], 
the behaviour of N{r) at small r differs for the cases of fixed and running coupling. 
For small r < l/Qs{Y), forms of the type (rQs)^'^ \n{CrQs) [33] describe the 
solutions at sufficiently high rapidity, where 7, defined in a y-independent fitting 
region, is ~ 0.65 for the fixed coupling constant but 7 ~ 0.85 for running coupling. 
These values are for the limit Y 00. 

4. Arguments in [32,33] suggest a lower limit to the scaling window rQs{Y) ~ 
AQCD/Qs{y) below which N{r) returns to the perturbative double-leading-loga- 
rithmic expression. Remarkably, the scaling forms (fT^ proposed in [33, 45] give 
good fits to the solutions of BK even outside the scaling window, for r < A/Ql{Y), 
A ~ 0.2 GeV. Hence, it is not possible to establish numerically the limit of 
the scaling region as a deviation from scaling. However, the double-leading-log 
approximation provides an equally good description of the numerical solution in 
the r-region below the scaling window. 

5. For fixed coupling, the scale invariance of the kernel preserves any A-dependence 
of the initial condition during BK evolution. For running coupling and for very 
large energies and nuclear sizes, we have re-derived and checked numerically 
Equation the A-dependence decreases with increasing rapidity like l/y/Y 
[34]. 

The above results have been established by evolving over many orders of magnitude 
in energy. Thus, any phenomenological application of these findings has to assume that 
initial conditions can be fixed at (and perturbatively evolved from) a sufficiently small 
energy scale for the non-linear evolution to be effective in an experimentally accessible 
regime. Moreover, phenomenology based on the BK equation will face at least some 
of the problems known from applications of BFKL such as the question of whether 
and how to implement kinematical cuts for gluon emission. Despite these caveats, it 
is interesting to compare the numerical results found here to the general trends in 
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the data. A comparison of saturation-inspired parameterizations with data on lepton- 
proton, lepton-nucleus and nuclear coUisions at high energies suggests a saturation 
scale oc A° exp{DY) with D ~ 0.29 [88] and a ~ 4/9 > 1/3 [57] (for related 
phenomenological studies, see [55,56]). 

Our results allow to discuss to which extent existing data, showing geometric scal- 
ing, differ from the asymptotic BK scaling behaviour. In particular, the strong A- 
dependence of the saturation scale seen in the data indicates, at variance with the 
result from the BK scaling solution with running coupling, that the properties of the 
initial nuclear condition have not yet been washed out by non-linear small-x evolu- 
tion. The kinematic range of the lepton-nucleus data studied in [56, 57] is too small to 
test this evolution. Also, the exponential F-dependence of the saturation scale with 
D ~ 0.3 seen in the data can appear naturally from BK evolution of reasonable initial 
conditions over some units in rapidity in the running coupling case. But this value of D 
is not a property of the asymptotic solution for running coupling. For fixed coupling, 
it can only be obtained with unrealistically small values of the coupling constant. 

None of these facts contradicts non-hnear BK evolution - they simply illustrate 
that the evolution observed in experimental data has not yet reached its asymptotic 
behaviour. To further advance our understanding of saturation effects in QCD dynam- 
ics at high energies, both theoretical and experimental studies are required. In the 
context of the BK equation, this requires the study of solutions under more realistic 
conditions. In particular, the impact parameter dependence may have a significant 
effect on the A-dependence of the saturation scale, a point which we plan to study in 
the future. On the experimental side, the forward rapidity measurements at the BNL 
Relativistic Heavy Ion CoUider RHIC give access to a kinematic window interesting 
for small-x evolution studies. These studies are at the very beginning. Also, in the 
near future measurements at the CERN Large Hadron Collider LHC will provide more 
stringent tests of small-x evolution, extending the kinematic reach by at least three 
orders of magnitude further down in the momentum fraction x. 
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